Scroll to
  1. Simulation of N particle Interaction via Lennard Jones potential using python
  2. Abstract/Aim
  3. 1Introduction
  4. 1.1Lennard Jones Potential
  5. 1.2N particle Interaction
  6. 2Method
  7. 2.1Algorithm
  8. 2.1.1Predictor Corrector method
  9. 2.1.2Runge Kutte method
  10. 2.2Interaction1
  11. 2.2.1All pair method
  12. 2.2.2Neighbour List method
  13. 2.3Time saving techniques
  14. 2.4Setting the parameters
  15. 2.5Law of Conservation Of Energy as test of Reliabilty
  16. 2.6Coding Essentials
  17. 2.6.1Object Oriented Programming
  18. 2.7Code
  19. 3Result
  20. 3.1Initial T=79 K , ρ=1.3954  g/cc\rho=1.3954\;g/ccρ=1.3954g/cc<math xmlns="http://www.w3.org/1998/Math/MathML"><mi>&#x3C1;</mi><mo>=</mo><mn>1</mn><mo>.</mo><mn>3954</mn><mo>&#xA0;</mo><mi>g</mi><mo>/</mo><mi>c</mi><mi>c</mi></math>\rho=1.3954\;g/cc
  21. 3.2Initial T=163 K , ρ=1.3954  g/cc\rho=1.3954\;g/ccρ=1.3954g/cc<math xmlns="http://www.w3.org/1998/Math/MathML"><mi>&#x3C1;</mi><mo>=</mo><mn>1</mn><mo>.</mo><mn>3954</mn><mo>&#xA0;</mo><mi>g</mi><mo>/</mo><mi>c</mi><mi>c</mi></math>\rho=1.3954\;g/cc
  22. 3.3Initial T=202 K, ρ=1.3954  g/cc\rho=1.3954\;g/ccρ=1.3954g/cc<math xmlns="http://www.w3.org/1998/Math/MathML"><mi>&#x3C1;</mi><mo>=</mo><mn>1</mn><mo>.</mo><mn>3954</mn><mo>&#xA0;</mo><mi>g</mi><mo>/</mo><mi>c</mi><mi>c</mi></math>\rho=1.3954\;g/cc
  23. 3.4Initial T=98 K, ρ=0.697  g/cc\rho=0.697\;g/ccρ=0.697g/cc<math xmlns="http://www.w3.org/1998/Math/MathML"><mi>&#x3C1;</mi><mo>=</mo><mn>0</mn><mo>.</mo><mn>697</mn><mo>&#xA0;</mo><mi>g</mi><mo>/</mo><mi>c</mi><mi>c</mi></math>\rho=0.697\;g/cc​
  24. 3.5Initial T=168 , ρ=0.697  g/cc\rho=0.697\;g/ccρ=0.697g/cc<math xmlns="http://www.w3.org/1998/Math/MathML"><mi>&#x3C1;</mi><mo>=</mo><mn>0</mn><mo>.</mo><mn>697</mn><mo>&#xA0;</mo><mi>g</mi><mo>/</mo><mi>c</mi><mi>c</mi></math>\rho=0.697\;g/cc​
  25. 3.6Initial T=214 K , ρ=0.697  g/cc\rho=0.697\;g/ccρ=0.697g/cc<math xmlns="http://www.w3.org/1998/Math/MathML"><mi>&#x3C1;</mi><mo>=</mo><mn>0</mn><mo>.</mo><mn>697</mn><mo>&#xA0;</mo><mi>g</mi><mo>/</mo><mi>c</mi><mi>c</mi></math>\rho=0.697\;g/cc​
0%

Simulation of N particle Interaction via Lennard Jones potential using python

Harsimran Singh

B.Sc-(Honours)-Physics-(III), Panjab University , Chandigarh

Abstract/Aim

Simulation of N particle interaction has been key for many advances in fields like Solid State Physics, Molecular Dynamics, Astrophysics and much more. As many times, systems under observation do not contain one atom or molecule but many of them arranged in random way where each configuartion has its own probability of occurence. In all such system one has to consider N body interaction inevitably. Due to such profound occurence, N body interaction simulation is strong tool which can be used in many fields like plasma physics. One use of N body interaction simulation had been in study of behaviour of gases and estimation of properties of atoms from behaviour of bulk. One of its example has been simulation of behaviour of argon which has been explained quite well by lennard jones potential. Due to its short range nature, It is the potential for which techniques of Molecular dynamics simulation works really well and it's computational cost is less. In this project, Simulation for Argon atoms interacting via l-J potential has been done using python3.

Introduction

Lennard Jones Potential

Lennard jones potential is short range potential which is attractive over large distancces but repulsive at very short distance. Form of this potential is

V=4ε((σr)12(σr)6)

In this , minimum valuue of potential is at r=21/6σ and Vlj=ε

lennard Jones potential is the one of the most used potential while modeling for basic properties of neutral atoms and molecules, and both attractive and repulsive parts which account for van der waal forces and interatomic repulsion respectively . One of the most important advantage of using lj potential is because it is short range potential and dies out rapidly, which in turn opens the gateway for one simulation technique which makes computation times much lower than that of any other long range potential.

potential.png
     Lennard jones potential as function of r

    N particle Interaction

    In system containing more than one or two atoms or molecules, one has to consider interaction of each molecule with N-1 other atoms. From perspective of quantum mechanics, this is translated by Schroedinger wave equation.

    iψt=Hψ

    here hamiltonion include potential due to all molecules and psi represent state of system. In classical mechanics, equation governing time evolution of system can be written as

    d2ridt2=1miijjfij

    For simulation and computation we will be using classical equation and lj potential as governing equation of motion.

    F=Vlj

    Method

    Main aim of doing simulation is to compute position r and velocity v of each atom for later time t. this is done by computing position and velocity after each time dt and then repeating this process untill we reach desired state or time. Main part of program is to calculation of r(t+h) of each atom from previous configuation. Here we used two algorithm's. However only one of them has been used mostly used while other was used for confirmation of stability and uniqueness of solution.

    Algorithm

    Predictor Corrector method

    Predictor corrector method is also known as modified euler method. It has been used excessively in domain of molecular dynamics because of its simplicity, reliability and lesser computation cost. ​This method had been used for simulating Argon behaviour (Rahman, et al, 1964 ). this algorithm is as:

    Assuming we know position of an atom at time t=t-h , which we write as r n-1 , and position , velocity, accelaration at t=t, which is r n , v n, a n, now first we predict r n+1 and then refine it later on as

    ren+1=rn1+2tvn

    ​from estimated position , we compute estimated accelaration ae n+1. After this, we calculate corrected value of r and v as

    vn+1=vn+12t(an+aen+1)
    rn+1=rn+12t(vn+vn+1)

    a n+1 can now be computed from calculated r n+1. and process goes on.

    Runge Kutte method

    RK method is quite known method, But due to more computational cost it is only used here as tester for stability and uniquesness of solution provided by predictor corrector method.

    k0=t2f(rn)
    k1=t2f(rn+12tvn+18k0)
    k2=t2f(rn+tvn+12k1)
    rn+1=rn+tvn+16(k0+2k1)
    vn+1=vn+16t(k0+4k1+k2)

    Interaction1

    All pair method

    It is the simplest to implement, but extremely inefficient when the interaction range r c is small compared with the linear size of the simulation region. All pairs of atoms must be examined because, owing to the continual rearrangement that characterizes the fluid state, it is not known in advance which atoms actually interact. Although testing whether atoms are separated by less than r c is only a part of the overall interaction computation, the fact that the amount of computation needed grows as O(N 2) rules out the method for all but the smallest values of N . Techniques for reducing this growth rate to a more acceptable O(N ) level, often used in tandem, will be discussed here. A schematic summary of the methods appears in Figure 2.

    Screenshot from 2021-04-30 00-00-40.png
      All pair method

      Neighbour List method

      In case of N particle interaction, computatiions are required for N(N-1)/2 pairs which results in computation varying roughly as N2. For decreasing this dependence , In case of short range potential , assign neighbour list to all atoms , by which only those present in neighbours are considered for force calculation and this list can be updated after 10-20 cycles which makes computation time vary roughly as N and not as N2.

      Screenshot from 2021-04-29 22-44-48.png
        figure depicting neighbor list method

        Time saving techniques

        To reduce time of computation which has been major drawbac of python3 , different compiler based python was used named pypy. It proved to make code run 6 times faster than on normal python, but was not capable of plotting and animation and thus all the data was once computed and saved and then was plottted and animate dusing python3.

        Setting the parameters

        for Argon, mass of an atom is 39.95x1.67x10-27 Kg. In Lennard Jones Potential, we have σ=3.4  Angstorm(​Rahman, et al, 1964),​(Verlet, et al, 1967​).

        In order to make compuatation easier following units have been adopted:

        length-3.4 Amgstorm ( σ)

        Time- 10-14 sec​

        ε=KbTa(Ta=120K)

        from now on we will call our length unit X and time unit C (not to be confused by speed of light).

        In these units, Boltzmann Constant is

        Kboltzmann=1.064 x 10-31 Kg X2/C2

        Force equation in pair is

        F=24  ε  1r2(2r121r6)((xx0)i+(yy0)j+(zz0)k)

        Thus This equation becomes:

        d2rdt2=24  εM  1r2(2r121r6)((xx0)i+(yy0)j+(zz0)k)

        In our present units, for Argon, This eqquation becomes

        d2rdt2=0.0004593  1r2(2r121r6)((xx0)i+(yy0)j+(zz0)k)

        Law of Conservation Of Energy as test of Reliabilty

        In Simulation, One is not aware of final output and depends entirely on simulation itself for validation. In such cases, Whether our Algorithm has errors under control or not is hard to determine in exact way because one do not know theoretical values of such system. In order to accomplish this, Method of Law of Conservation of Energy is used. System Initial total energy is compared with system final Energy to get the Idea about Uniqueness and stability of Code.

        Coding Essentials

        Code has been made using Object Oriented Programming Paradigm. It is due to the fact that we need to take care of properties of particles and which are huge in number , thus makig it quite redundant to implement in procedural programming. As, Python is object oriented language along with its ability to loop over objects, It became primary choice. Some basics are discussed here.

        Object Oriented Programming

        Object oriented programming allows us to make objects with specific attributes and methods, which we used to create particles. As our particles have properties like position, velocity, accelaration associated making a class particle was inevitabe part. Here we discuss brief method of making class in python3.

        1. Class Particle(): ///// Starting Class Definition
        2. def __init__(self,position,velocity):

        Here position and velocity are additional Arguments, which one hould enter while creating an object.

        1. class particle():
        2. def __init__(self, position,velocity):
        3. self.position=position /// Defining attributes of object
        4. self.velocity=velocity
        5. self.accelaration=(0,0,0)

        Here self.X is attribute with name X which belongs to the object , in our case particle. As it is evident that every Object attribute need not to be in definition of __init__().

        Now we will se how to add function in our class and objects. Here Two distinction occurs. In code two kind of Methods(Functions in Class are known as methods).

        1. Methods which belong to objects. That is every object have access to them. They are defined as normal functions.

        1. def method_name(self,*Args): // self is our object. , *Args are other user defined arguments
        2. ///Using This method
        3. object.method_name(*Args)

        2. Methods which belong to class only and cannot be accessed by objects are Classmethods. They are also defined in same way except we have to specify that it is class method

        1. @classmethod
        2. def classmethod_name(cls, *Args): \\\ cls is argument used to denote class, *Args are other user defined arguments.
        3. \\\Using class Method. In our case class name is Particle.
        4. Particle.classmethod_name(*Args)

        For More Info about Python classes and objects. See-(Learn Python-The Hard way)

        Code

        we first define class particle which was used to create particle objects with attributes required in computation for each particle ,object and class methods.

        Defining Class

        1. ### Defining class particle for creating particle objects
        2. class particle():
        3. ### List to store handle of each particle created
        4. objects = []
        5. # time variable for keeping track of time of simulation
        6. time=0
        7. # defining initiating function of class, which creates object
        8. def __init__(self, position, velocity):
        9. #position, velocity , accelaration and neighbour list attribute of particle
        10. self.position = position
        11. self.velocity = velocity
        12. self.accelaration = (0, 0, 0)
        13. self.neighbours = []
        14. # temporary position and velocity variable used in algorithm
        15. self.position1=position
        16. self.velocity1=velocity
        17. self.velocity2=velocity
        18. self.velocity3=velocity
        19. self.accelaration1=(0,0,0)
        20. self.da=(0,0,0)
        21. # adding object handle in objects list
        22. particle.objects.append(self)
        23. # method for updating position and velocity by 1 step
        24. @classmethod
        25. def update(cls):
        26. # step size
        27. t=0.1
        28. # half of step size required for estimation
        29. th=t/2
        30. # updating neighbours
        31. if particle.time%100==0:
        32. particle.update_neighbours()
        33. # updating accelaration before updating coordinates
        34. particle.update_ac()
        35. # calculating position and velocity
        36. for obj in cls.objects:
        37. # position at n-1 step
        38. xn1, yn1, zn1 = obj.position1
        39. # velocity
        40. vx,vy,vz=obj.velocity
        41. # creating box by condition
        42. if xn1>=2.7:
        43. vx,vy,vz=(-abs(vx),vy,vz)
        44. if yn1>=2.7:
        45. vx,vy,vz=(vx,-abs(vy),vz)
        46. if zn1>=2.7:
        47. vx,vy,vz=(vx,vy,-abs(vz))
        48. if xn1<=-2.7:
        49. vx,vy,vz=(abs(vx),vy,vz)
        50. if yn1<=-2.7:
        51. vx,vy,vz=(vx,abs(vy),vz)
        52. if zn1<=-2.7:
        53. vx,vy,vz=(vx,vy,abs(vz))
        54. obj.velocity=(vx,vy,vz)
        55. # storing current velocity and accelaration and position in temporary variable
        56. obj.velocity1=obj.velocity
        57. obj.accelaration1=obj.accelaration
        58. obj.position1=obj.position
        59. # updating poition for estimated value at n+1 step
        60. obj.position = (xn1 + 2*t * vx, yn1 + 2*t * vy, zn1 + 2*t * vz)
        61. # updating accelaration for estimated location
        62. particle.update_ac()
        63. # correcting estimate
        64. for obj in particle.objects:
        65. vxn,vyn,vzn=obj.velocity
        66. obj.velocity=(vxn+(obj.accelaration[0]+obj.accelaration1[0])*th,vyn+(obj.accelaration[1]+obj.accelaration1[1])*th,vzn+(obj.accelaration[2]+obj.accelaration1[2])*th)
        67. obj.position=(obj.position1[0]+(vxn+obj.velocity[0])*th,obj.position1[1]+(vyn+obj.velocity[1])*th,obj.position1[2]+(vzn+obj.velocity[2])*th)
        68. # incrementing time
        69. particle.time+=1
        70. # defining force function
        71. def force(self):
        72. ep = 0.0004593 # Well depth kb_m*T , T= 120K
        73. d = 1 # sigma
        74. f = (0, 0, 0)
        75. x0, y0, z0 = self.position
        76. for obj2 in self.neighbours:
        77. x, y, z = obj2.position
        78. r=((x-x0)**2+(y-y0)**2 +(z-z0)**2)
        79. if r<7.84:
        80. p = ep * (2 * ((1)/(r**7)) - 1 * ((1)/(r**4)))
        81. f = (f[0] -p * (x - x0), f[1] - p * (y - y0), f[2] - p * (z - z0))
        82. # changing accelaration so as to not calculate force for future config.
        83. self.accelaration=f
        84. # creating update neighbour function
        85. @classmethod
        86. def update_neighbours(cls):
        87. R_critical=3 # 2.5*Sigma
        88. for obj in cls.objects:
        89. obj.neighbours=[]
        90. x,y,z=obj.position
        91. for obj2 in cls.objects:
        92. x1,y1,z1=obj2.position
        93. r=((x-x1)**2+(y-y1)**2+(z-z1)**2)**0.5
        94. if r<R_critical and r>0:
        95. obj.neighbours.append(obj2)
        96. @classmethod
        97. def update_neighbours1(cls):
        98. R_critical=2.8 # 2.5*Sigma
        99. for obj in cls.objects:
        100. obj.neighbours=[]
        101. x,y,z=obj.position
        102. for obj2 in cls.objects:
        103. x1,y1,z1=obj2.position
        104. r=((x-x1)**2+(y-y1)**2+(z-z1)**2)**0.5
        105. if r<R_critical and r>0:
        106. obj.neighbours.append(obj2)
        107. # calculating energy
        108. @classmethod
        109. def energy(cls):
        110. e = 0
        111. for obj in cls.objects:
        112. vx, vy, vz = obj.velocity
        113. e = e + vx**2 + vy**2 + vz**2
        114. return e
        115. # creating function to update acceleration
        116. @classmethod
        117. def update_ac(cls):
        118. for obj in cls.objects:
        119. obj.force()
        120. x_list = []
        121. y_list = []
        122. z_list = []
        123. # function to get position of all particles
        124. @classmethod
        125. def getcoordinate(cls):
        126. cls.x_list = []
        127. cls.y_list = []
        128. cls.z_list = []
        129. for obj in cls.objects:
        130. x, y, z = obj.position
        131. cls.x_list.append(x)
        132. cls.y_list.append(y)
        133. cls.z_list.append(z)
        134. return cls.x_list, cls.y_list, cls.z_list
        135. # function to get velocities of all particles
        136. @classmethod
        137. def getvelocity(cls):
        138. cls.x_list = []
        139. cls.y_list = []
        140. cls.z_list = []
        141. for obj in cls.objects:
        142. x, y, z = obj.velocity
        143. cls.x_list.append(x)
        144. cls.y_list.append(y)
        145. cls.z_list.append(z)
        146. return cls.x_list, cls.y_list, cls.z_list
        147. # functions for generating velocity distribution data initial and final
        148. @classmethod
        149. def velocity_distribution(cls):
        150. v=[]
        151. for obj in cls.objects:
        152. v.append((obj.velocity[0]**2+obj.velocity[1]**2+obj.velocity[2]**2)**0.5)
        153. np.save('velocity_data.npy',v,allow_pickle=True)
        154. @classmethod
        155. def velocity_distribution1(cls):
        156. v1=[]
        157. for obj in cls.objects:
        158. v1.append((obj.velocity[0]**2+obj.velocity[1]**2+obj.velocity[2]**2)**0.5)
        159. np.save('velocity_data1.npy',v1,allow_pickle=True)
        160. # function for calculating potential energy due to neighbours
        161. @classmethod
        162. def potential(cls):
        163. potent=0
        164. for obj in cls.objects:
        165. x0,y0,z0= obj.position
        166. for obj2 in obj.neighbours:
        167. x, y, z = obj2.position
        168. r=((x-x0)**2+(y-y0)**2 +(z-z0)**2)**(0.5)
        169. potent+=480*(1/(r**12) - 1/(r**6))
        170. potent=potent/2
        171. return potent
        172. # function to calculate potential energy due to all particles
        173. @classmethod
        174. def potential1(cls):
        175. potent=0
        176. for obj in cls.objects:
        177. x0,y0,z0= obj.position
        178. for obj2 in cls.objects:
        179. x, y, z = obj2.position
        180. r=((x-x0)**2+(y-y0)**2 +(z-z0)**2)**(0.5)
        181. if r>0:
        182. potent+=480*(1/(r**12) - 1/(r**6))
        183. return potent/2
        defining class

        In order to create particles with given position and random velocity, create function was created as follows:

        1. def create():
        2. lists = []
        3. v_da1=[]
        4. x_da1=[]
        5. for i in range(125):
        6. x = 0.001*np.random.randint(-1000,1000)
        7. y = 0.001*np.random.randint(-1000,1000)
        8. z = 0.001*np.random.randint(-1000,1000)
        9. vx = 0.0001*rd.randint(-100,100)*0.867
        10. vy = 0.0001*rd.randint(-100,100)*0.867
        11. vz = 0.0001*rd.randint(-100,100)*0.867
        12. x_da1.append((x,y,z))
        13. v_da1.append((vx,vy,vz))
        14. lists.append(particle((x, y, z), (vx, vy, vz)))
        15. np.save("v_da1.npy",v_da1)
        16. np.save("x_da1.npy",x_da1)
        Create function

        In order to calculate local density, whole volume was divided into 64 parts. and number of particle inside that box were counted and used in calculating standard deviation of number of particles in each box. and this deviation was calculated for every time.

        1. import math
        2. density=np.zeros((4,4,4))
        3. density_data=[]
        4. for data in position:
        5. density=np.zeros((4,4,4))
        6. x,y,z=data
        7. for number in range(len(x)):
        8. density[1+math.ceil(x[number]/7.51),1+math.ceil(y[number]/7.51),1+math.ceil(z[number]/7.51)]+=1
        9. density_data.append(np.std(density)/np.mean(density))
        10. np.save('density.npy',density_data,allow_pickle=True)
        Denisty Deviation

        At the end we have animation and plotting code

        1. import os
        2. import numpy as np
        3. print('\n Plotting and animation Begins\n\n')
        4. os.chdir('')
        5. a=np.load('data_set.npy',allow_pickle=True)
        6. v1=np.load('velocity_data1.npy',allow_pickle=True)
        7. v=np.load('velocity_data.npy',allow_pickle=True)
        8. T_data=np.load('vel_data.npy',allow_pickle=True)
        9. step=0.0005
        10. final=0.02
        11. mean_v=np.mean(v)
        12. mean_v1=np.mean(v1)
        13. meansv=np.mean(v**2)**0.5
        14. meansv1=np.mean(v1**2)**0.5
        15. M=39.95*1.67*10**(-27)
        16. Kb=1.064*10**(-32)
        17. #Temperature
        18. Ti=(M/(3*Kb))*(np.mean(v1**2))
        19. Tf=(M/(3*Kb))*(np.mean(v**2))
        20. vr=np.arange(0,final,step)
        21. Maxw=125*(4*np.pi*vr**2*step)*((M/(2*np.pi*Kb*Tf))**(1.5))*np.exp(-(M/(2*Kb*Tf))*vr**2)
        22. print('Mean v:',np.mean(v))
        23. import matplotlib.pyplot as plt
        24. plt.plot(T_data)
        25. plt.title('Variation Of temperature with Time')
        26. plt.ylabel('Temparature')
        27. plt.ylim(60,1.1*np.max(T_data))
        28. plt.figure()
        29. plt.subplot(121)
        30. a1,b1,c1=plt.hist(v1,bins=np.arange(0,final,step))
        31. ind=np.argmax(a1)
        32. plt.title("Initial Velocity Distribution")
        33. plt.subplot(122)
        34. a2,b2,c2=plt.hist(v,bins=np.arange(0,final+step,step))
        35. ind2=np.argmax(a2)
        36. if np.max(a1)>np.max(a2):
        37. limit=np.max(a1)
        38. else:
        39. limit=np.max(a2)
        40. plt.title("Final Velocity Distribution")
        41. plt.text(0.012,1.4*limit-4,"Temperature: {:.2f}".format(Tf), bbox=dict(boxstyle='round',facecolor='white'))
        42. plt.text(0.012,1.4*limit-5,"Mean v : {:.4f}".format(mean_v)+"X/c", bbox=dict(boxstyle='round',facecolor='white'))
        43. plt.text(0.012,1.4*limit-6,"RMS v : {:.4f}".format(meansv)+"X/c", bbox=dict(boxstyle='round',facecolor='white'))
        44. plt.text(0.012,1.4*limit-7,"MP v : {:.4f}-{:.4f}".format(b2[ind2],b2[ind2+1])+"X/c", bbox=dict(boxstyle='round',facecolor='white'))
        45. plt.plot(vr,Maxw)
        46. plt.legend(["Maxwell Distribution","Observed Distribution"])
        47. plt.subplot(121)
        48. plt.ylim(0,1.4*limit)
        49. plt.text(0.012,1.4*limit-4,"Temperature: {:.2f}".format(Ti), bbox=dict(boxstyle='round',facecolor='white'))
        50. plt.text(0.012,1.4*limit-5,"Mean v : {:.4f}".format(mean_v1)+"X/c", bbox=dict(boxstyle='round',facecolor='white'))
        51. plt.text(0.012,1.4*limit-6,"RMS v : {:.4f}".format(meansv1)+"X/c", bbox=dict(boxstyle='round',facecolor='white'))
        52. plt.text(0.012,1.4*limit-7,"MP v : {:.4f}-{:.4f}".format(b1[ind],b1[ind+1])+"X/c", bbox=dict(boxstyle='round',facecolor='white'))
        53. plt.subplot(122)
        54. plt.ylim(0,1.4*limit)
        55. plt.show()
        56. print('###################')
        57. density=np.load('density.npy',allow_pickle=True)
        58. plt.plot(density)
        59. plt.title('Variation of Deviation in Local Density with time')
        60. plt.ylim(0.1,1.1*np.max(density))
        61. plt.show()
        62. print('wait a second')
        63. def update_coordinate(i, scatter):
        64. x, y, z = a[i]
        65. ax.clear()
        66. ax.set_xlim3d([-3,3])
        67. ax.set_ylim3d([-3,3])
        68. ax.set_zlim3d([-3,3])
        69. ax.set_title('Interaction Animation')
        70. scatter[0]=ax.scatter3D(x,y,z,color='green',marker=".")
        71. import mpl_toolkits.mplot3d.axes3d as mp
        72. import matplotlib.animation as animation
        73. Writer = animation.writers['ffmpeg']
        74. writer = Writer(fps=10, metadata=dict(artist='Me'), bitrate=2800)
        75. fig=plt.figure()
        76. ax = plt.axes(projection ="3d")
        77. ax.set_xlim3d([-3,3])
        78. ax.set_ylim3d([-3,3])
        79. ax.set_zlim3d([-3,3])
        80. ax.set_title("Scatter Animation")
        81. x,y,z=a[0]
        82. scatter=[ax.scatter3D(x,y,z,color='green')]
        83. anim=animation.FuncAnimation(fig,update_coordinate,len(a),fargs=(scatter,),interval=20)
        84. anim.save("LJ_Potential.mp4",writer=writer)
        85. plt.show()
        Animation and plot code

        one script file was written to run all codes on different platform as required.

        1. pypy3 main.py
        2. pypy3 density_distribution.py
        3. python3 animate.py
        Script file for terminal

        Result

        Simulation was carried out for two different Densities of Argon p=1.3954 g/cc and p=0.6977 g/cc. foe three different Temperatures each.

        Initial T=79 K , ρ=1.3954  g/cc

        Initial Temperature=79 K
         Step Size=0.1
         Integration Time: 50,000 x 0.1 
        Time taken for execution : 205.84 sec 
         Initial average Kinetic Energy=118.67 Kb(Boltzmann Constant)
        Initial average potential energy of an atom = -221 Kb
        Final Kinetic Energy= 133.88 Kb
         Final average potential energy of an atom = -235.44 Kb
         Total error in Temperature = 0.5153 K
        Total Error in energy in Kb = 0.773 
        Animation 
          T80pcD_1.png
          Velocity Distribution
            T80pcT.png
            Variation Of Temperature with time
              T80pcL.png
              Variation of local density deviation with time
                Results of simulation

                Initial T=163 K , ρ=1.3954  g/cc

                Initial Temperature=163.3 K
                 Step Size=0.1
                 Integration Time: 50,000 x 0.1 
                Time taken for execution : 50.95 sec 
                 Initial average Kinetic Energy=244.95 Kb(Boltzmann Constant)
                Initial average potential energy of an atom = -221 Kb
                Final Kinetic Energy= 136.82 Kb
                 Final average potential energy of an atom = -106.76 Kb
                 Total error in Temperature = 4.073 K
                Total Error in energy in Kb = 6.1105 
                Animation
                  T163pcD.png
                  Velocity Distribution
                    T163pcT.png
                    Variation of Temperature with Time
                      T163pcL.png
                      Variation of Local density deviation with time

                        Initial T=202 K, ρ=1.3954  g/cc

                        Initial Temperature=202.3 K
                         Step Size=0.1
                         Integration Time: 50,000 x 0.1 
                        Time taken for execution : 29.19 sec 
                         Initial average Kinetic Energy=303.5 Kb(Boltzmann Constant)
                        Initial average potential energy of an atom = -221 Kb
                        Final Kinetic Energy= 140.91 Kb
                         Final average potential energy of an atom = -49.62 Kb
                         Total error in Temperature = 5.86 K
                        Total Error in energy in Kb = 8.739
                        Animation
                          T202pcD.png
                          Velocity Distribution
                            T202pcT.png
                            Variation of Temperature with Time
                              T202pcL.png
                              Variation of Local density deviation with time

                                Initial T=98 K, ρ=0.697  g/cc

                                Initial Temperature=97.6 K
                                 Step Size=0.1
                                 Integration Time: 50,000 x 0.1 
                                Time taken for execution : 19.16 sec 
                                 Initial average Kinetic Energy=146.41 Kb(Boltzmann Constant)
                                Initial average potential energy of an atom = -60.37 Kb
                                Final Kinetic Energy= 117.93 Kb
                                 Final average potential energy of an atom = -26.16 Kb
                                 Total error in Temperature = 3.81 K
                                Total Error in energy in Kb = 5.72
                                Animation
                                  T98ptD.png
                                  Velocity Distribution
                                    T98ptT.png
                                    Variation of Temperature with Time
                                      T98ptL.png
                                      Variation of Local density deviation with time

                                        Initial T=168 , ρ=0.697  g/cc

                                        Initial Temperature=168.3 K
                                         Step Size=0.1
                                         Integration Time: 50,000 x 0.1 
                                        Time taken for execution : 14.60 sec 
                                         Initial average Kinetic Energy=252.8 Kb(Boltzmann Constant)
                                        Initial average potential energy of an atom = -60.37 Kb
                                        Final Kinetic Energy= 211.75 Kb
                                         Final average potential energy of an atom = -13.41 Kb
                                         Total error in Temperature = 3.93 K
                                        Total Error in energy in Kb = 5.91
                                          T168ptD.png
                                          Velocity Distribution
                                            T168ptT.png
                                            Variation of Temperature with Time
                                              T168ptL.png
                                              Variation of Local density deviation with time

                                                Initial T=214 K , ρ=0.697  g/cc

                                                Initial Temperature=214.29 K
                                                 Step Size=0.1
                                                 Integration Time: 50,000 x 0.1 
                                                Time taken for execution : 13.34 sec 
                                                 Initial average Kinetic Energy=321.44 Kb(Boltzmann Constant)
                                                Initial average potential energy of an atom = -60.37 Kb
                                                Final Kinetic Energy= 280.70 Kb
                                                 Final average potential energy of an atom = -13.84 Kb
                                                 Total error in Temperature = 3.87 K
                                                Total Error in energy in Kb = 5.80
                                                Animation
                                                  T214ptD.png
                                                  Velocity Distribution
                                                    T214ptT.png
                                                    Variation of Temperature with Time
                                                      T214ptL.png
                                                      Variation of Local density deviation with time

                                                        References

                                                        • Rahman, A. (1964) “Correlations in the Motion of Atoms in Liquid Argon.”

                                                        • Verlet, Loup (1967) “Computer ‘Experiments’ on Classical Fluids. I. Thermodynamical Properties of Lennard-Jones Molecules.”

                                                        Source

                                                        • Fig 2: D.C. Rapport, Cambridge University Press, The Art of molecular Dynamics Simulation
                                                        • Fig 3: The Art of Molecular Dynamics by D.C. Rapport
                                                        Manage ContentSee Navigation
                                                        • 0 Simulation of N particle Interaction via Lennard Jones potential using python
                                                        • 0 Abstract/Aim
                                                        • 0 1 Introduction
                                                        • 0 1.1 Lennard Jones Potential
                                                        • 0 1.2 N particle Interaction
                                                        • 0 2 Method
                                                        • 0 2.1 Algorithm
                                                        • 0 2.1.1 Predictor Corrector method
                                                        • 0 2.1.2 Runge Kutte method
                                                        • 0 2.2 Interaction1
                                                        • 0 2.2.1 All pair method
                                                        • 0 2.2.2 Neighbour List method
                                                        • 0 2.3 Time saving techniques
                                                        • 0 2.4 Setting the parameters
                                                        • 0 2.5 Law of Conservation Of Energy as test of Reliabilty
                                                        • 0 2.6 Coding Essentials
                                                        • 0 2.6.1 Object Oriented Programming
                                                        • 0 2.7 Code
                                                        • 0 3 Result
                                                        • 0 3.1 Initial T=79 K , ρ=1.3954  g/cc\rho=1.3954\;g/cc
                                                        • 0 3.2 Initial T=163 K , ρ=1.3954  g/cc\rho=1.3954\;g/cc
                                                        • 0 3.3 Initial T=202 K, ρ=1.3954  g/cc\rho=1.3954\;g/cc
                                                        • 0 3.4 Initial T=98 K, ρ=0.697  g/cc\rho=0.697\;g/cc
                                                        • 0 3.5 Initial T=168 , ρ=0.697  g/cc\rho=0.697\;g/cc
                                                        • 0 3.6 Initial T=214 K , ρ=0.697  g/cc\rho=0.697\;g/cc